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Abstract — Simulated landscapes have been used for decades 
to evaluate search strategies whose goal is to find the landscape 
location with maximum fitness. Applications include modeling 
the capacity of enzymes to catalyze reactions and the clinical 
effectiveness of medical treatments. Understanding properties of 
landscapes is important for understanding search difficulty. This 
paper presents a novel and transparent characterization of NK 
landscapes. 

We prove that NK landscapes can be represented by para- 
metric linear interaction models where model coefficients have 
meaningful interpretations. We derive the statistical properties 
of the model coefficients, providing insight into how the NK 
algorithm parses importance to main effects and interactions. An 
important insight derived from the linear model representation is 
that the rank of the linear model defined by the NK algorithm is 
correlated with the number of local optima, a strong determinant 
of landscape complexity and search difficulty. We show that 
the maximal rank for an NK landscape is achieved through 
epistatic interactions that form partially balanced incomplete 
block designs. Finally, an analytic expression representing the 
expected number of local optima on the landscape is derived, 
providing a way to quickly compute the expected number of 
local optima for very large landscapes. 

Index Terms — Balanced Incomplete Block Design, Orthant 
Probability, Walsh function 

I. Introduction 

SIMULATED landscapes have been used for several 
decades to evaluate search strategies whose goal is to find 
the landscape location with maximum fitness. Applications 
of simulated landscapes include modeling the capacity of 
enzymes to catalyze reactions or ligands to bind to proteins, 
and the clinical effectiveness of medical treatments [10], 0. 
Understanding properties of landscapes is important for un- 
derstanding search difficulty. 

NK landscapes are defined by a straightforward, tunable 
algorithm where N specifies the number of binary features 
or loci and K the degree of epistatic interactions among the 
loci [9]. NK landscapes are convenient because there are only 
two tunable parameters (N and K), and yet they provide a very 
rich set of landscapes. There is a large literature exploring 
properties of NK landscapes, though there are few analytic 
results, with J4| and flTl notable exceptions. 

This paper provides a novel perspective on NK landscapes 
and generalizations of NK landscapes, proving that these 
landscapes can be characterized by parametric linear models 
comprised of main effects and interaction effects where the 
model coefficients have meaningful interpretations. The NK 
algorithm induces a statistical distribution on the paramet- 
ric model coefficients. We derive the distribution of model 
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coefficients, showing how the NK algorithm, for K <C N, 
automatically assigns the largest expected magnitude to main 
effects, with the expected magnitude of interaction effects 
typically decreasing with increasing order of interaction. 

The linear model representation of the landscape suggests 
that model rank should provide a measure of landscape com- 
plexity. A simple method of assessing rank is provided, and we 
determine conditions on N and K sufficient for the existence 
of designs that achieve maximum rank. We then show that 
the expected number of local optima is proportional to an 
orthant probability, which can be calculated with reasonable 
speed and accuracy for very large landscapes. Interestingly and 
surprisingly, rank is both positively and negatively correlated 
with the number of local optima. For fixed N, it is well- 
known [3 1 that the number of local optima increases with K 
(i.e. a positive correlation with landscape rank). We show that 
when N and K are both fixed, there is a strong negative 
correlation between rank and number of local optima for 
classic NK landscapes. The statistical distribution of model 
effect coefficients provides an explanation for this counter 
intuitive phenomenon. 

In related work, a proposal to use linear models with main 
effects and interactions to construct landscapes was explored 
in fl4l . lfT31l and fl6l . These authors examined the effect of 
epistatic interactions on the properties of landscapes using 
metrics common in the experimental design literature, and they 
noted the equivalence between Walsh function decompositions 
and interaction models. An analysis of NK landscapes using 
Walsh functions was given in |6|. 

II. Generalized NK Landscapes and Interaction 
Models 

In this section we define the NK landscape and interaction 
models. The models are represented as linear models in matrix 
form, as this representation facilitates the study of model 
properties. For simplicity, the model is first discussed for K 
constant across loci. The model is then generalized to allow 
varying K across loci. 

A. NK landscapes 

A general landscape is defined by a triple (x, d, /) where 
X is a set of locations, d is a distance measure and / : x K 
is a "fitness" function. NK landscapes are a map from 
X = ({0, 1}) N to E where the fitness function / is built from 
N binary loci and epistatic interactions formed between each 
locus and K other loci where K can range from to N — 1. 

To describe NK landscapes more fully, for each i — 
1, . . . , 2 N let Xi € ({0, 1})^ represent i written in base 2 and 
represented as a (binary) vector of length N. For i = 1, . . . , N, 
let Wi denote a 2 K+1 x 1 vector of independent random 



2 



weights, where each component has mean fi/N and variance 
a 2 /N > and otherwise arbitrary probability distribution. 
Let fi(x) denote a function from ({0, 1}) N — > E 2 where 



E z 



denotes unit coordinate vectors in 



i.e. vectors 



of the form (0, 0,1,0,..., 0). The definition of the function 
fi(x) depends on the epistatic interactions to the ith locus. The 
function is defined explicitly below. 

For j = 1, ... , 2 N , define p 3 = J^^Li Wi T fi(xj). Note that 
Pj is the landscape fitness at location xj and that E[pj] = fi 
and Var(pj) = a 2 as each pj is comprised of a sum of N 
independent weights. The NK landscape is defined by the 
fitness-location pairs (pj,Xj) coupled with Hamming distance. 



B. Matrix representation of generalized NK landscapes 

We begin by providing an explicit description for the 
construction and representation of NK and generalized NK 
landscapes as linear models in matrix form. For i = 1, . . . , N, 
let Vi = {ii,z 2 ,...,iif i+ i} where i 2 , . . . , iKi+i) C 
{1, 2, . . . , N} and ii < i 2 < • • • < iRi+i with ij = i for some 
j. Vi denotes the ith interaction set, comprised of locus i and 
Ki loci that interact with locus i. Note that K isn't restricted 
to be constant across loci, and there are no restrictions on the 
number of times locus i can appear in the interaction sets. This 
represents a generalization of NK landscapes, with the classic 
NK landscape occurring as the special case when Ki = K 
and each locus appears in exactly K + 1 interaction sets. An 
additional generalization would be to not require locus i to be 
a member of Vi. The results in this article would still hold for 
this generalization. 

With a slight abuse of notation, for each i write fi(x) = 



_j ) where x h , . 



lK i+ 1 



are the elements 



of x corresponding to Vi. Thinking of {x ix , 
as a binary number, let e^x^, . . . ,%i K . +1 ) be the decimal 
representation of this number. Define the 1 x 2 Ki+1 vector 
fi(x n , . . . ,x iK +1 ) = (0, ...,0,1,0, ...,0) where the 1 oc- 



curs in column e^x^ . 



1 X lK i + 



1 ) + 1. Alternatively, if I 2 K i 



is the identity matrix of size 2 Ki+1 , then f^x^, . . . ,X{ K 
is the (ei(xi l: . . . ,Xi K +1 ) + l)th row of I 2 Ki+i. 

To write the generalized NK model in matrix form, consider 
the 1 x C vector f (x) = fx (x) | f 2 (x) | • • • | f N (x) where 
C = 2 Ki+1 an d I denotes column concatenation. Define 

the 2 N x C model matrix 



F 



( f(xi) \ 
f(x 2 ) 



\ f(X2N) / 



(1) 



With this definition, the generalized NK landscape is p = Fw 
where w = (wi T | • • • | wn t ) t denotes the Cxi vector of 
independent random fitnesses. 

Example: To illustrate the definition of and the matrix 
F, suppose N = 3, K % = 1 for i = 1, 2, 3 and Vi = {1, 2}, 
V 2 = {2, 3} and V 3 = {1, 3}. Then 



fi(xi,x 2 ) 



(1,0,0,0) if xx = 0, x 2 = 0; 

(0,1,0,0) if 2^=0, ir 2 = l; 

(0,0,1,0) if xi = l,z a = 0; 

(0,0,0,1) if xi = l,x 2 = 1. 



F 



The functions f 2 (x 2) X3) and f3(xi,a;3) are defined similarly. 
The matrix F is then 



/100010001000\ 
100001000100 
010000101000 
010000010100 
001010000010 
001001000001 
000100100010 

yoooioooioooi/ 

To our knowledge, NK landscapes have never been formal- 
ized using the matrix representation given here. The rank of 
F is a measure of the richness of the landscape, as it gives 
the dimension of the domain for the vector p. We will show 
that the rank of F is determined by N, {Ki}fL 1 and the 
structure of the interaction sets {Vi]f =1 , and that using rank 
as a measure of complexity provides refinement beyond using 
only N and {Ki\f =l . 

C. Interaction Model 

Here we define the general form of parametric interaction 
models, and in the next section relate them to NK landscapes. 
Statisticians have long employed interaction models to study 
the effects of multiple 'treatments' on an outcome, see for 
example lfT2l . Interaction models are straightforward to define 
and model parameters have meaningful interpretations. In 
the evolutionary computing literature, these models seem to 
have received little attention, with the exception of lfT4l , |fT5ll 
and lfT6l . 

To define a general interaction model, for k = 1, . . . , L, let 
Uk C {1,2,..., N} and define Xi = 2xi — 1. Mathematical 
properties of the interaction model and interpretation of model 
parameters are most easily obtained using the transformed 
values ii 6 {—1, 1} where as before Xi £ {0, 1}. 

The general form of an interaction model with L terms is 



k=l 



where the /3's are coefficients that can take any value in 
K and where we adopt the convention that when Uk = 0, 



n 



jeu k x j 



1. 



Example: Consider a model with N = 3 loci and U\ = 0, 

U 2 = 1,U 3 = 2,U 4 = 3,U 5 = {1,2},U 6 = {1,3}, U 7 = 
{2, 3}. The interaction model is q(x) = /3q> + f5\X\ + j3 2 x 2 + 
P3X3 + Pi 2 x 1 x 2 + /3 13 iix 3 + /3 23 i 2 x 3 . In this example, /3 2 
and 03 are main effects coefficients while Pi 2 ,(3i3 and (3 23 
are two-loci interaction coefficients. 
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The intercept (/3g) and coefficients of the main effects 
and interaction terms have meaningful interpretations. The 
intercept coefficient represents the average of the fitness values 
across the entire landscape. For the general interaction model, 
it is not difficult to show that the main effect 2/3^ represents 
the difference in fitness values when locus i is varied between 
Xi = 1 and Xi — — 1, averaged over the values of the 
other N — 1 loci. 4/J,j represents the difference of differences 
between fitness values for Xi = 1 and Xi = — 1 when Xj = 1 
and = — 1, averaged over all other Xk for k ^ In 
general, higher order interaction coefficients are interpreted 
as average differences between lower order interactions. For 
example, 8/J«fc represents the average difference in the two 
factor interaction between loci i and j when Xk = 1 and 
x fe = -1. 

III. Generalized NK Landscapes as Interaction 
Models 

Here we show that generalized NK landscapes can be 
expressed as linear interaction models, and that the interactions 
that are included in the model are completely determined by 
the interaction sets Vi, i = 1,...,N. We also derive the 
statistical properties of the interaction model coefficients and 
show how to construct classic NK landscapes that maximize 
the number of interaction terms. 

Let 2 Vi denote the power set for Vi, i.e. the set Vi and all 
it's subsets, including the empty set. Let T = UiLi 2 y \ and 
consider the interaction model 



p(x) = J2 II £ i- w 
c/gt jeu 

Evaluated at the 2^ transformed values xi,...,x 2 n, the 
model can be written in matrix notation as p = Fj3 where F 
is an appropriately defined 2 N x L matrix, L is the number of 
elements in T and j3 is an L x 1 vector of coefficients. The ran- 
dom vector j3 of main effects and interactions has distributional 
properties dependent on the probability distribution of the 
vector of weights w and the structure of the interaction sets Vi. 



The distributional properties of /3 are studied in section III-A 



Note that the model defined in Q always contains an intercept 
and all main effects terms. 



Example: Consider again the example in Section |II-B| 
where N = 3, Ki = 1 for i = 1,2,3 and V\ = 
{1,2}, V 2 = {2,3} and V 3 = {1,3}. Then T = 
({0},{1},{2},{3},{1,2},{2,3},{1,3}) and 



F = 



V i 



-1 
-1 
-1 
-1 
1 
1 
1 
1 



-1 
-1 
1 
1 

-1 
-1 
1 
1 



1 

-1 
-1 
1 

-1 

1 / 



I Pi \ 

02 
03 
012 
. @23 i 



The first column of F corresponds to the intercept /3g, columns 
two through four are for the main effects, and columns five 
through seven are for the two loci interactions. Column five, 
for example, corresponds to the interaction between loci 1 and 



2, and is obtained by taking the product of the elements of 
columns two and three, which correspond to the main effects 
for loci 1 and 2. 

The following proposition establishes that a generalized NK 
landscape defined by interaction sets Vy , . . . , Vm is equivalent 
to the interaction model given by pi, thereby clearly estab- 
lishing the nonzero interaction effects induced by the NK 
algorithm. The proof is given in the appendix. 

Proposition 1. Let F denote the model matrix for the gen- 
eralized NK landscape model defined in and F the 
model matrix for the interaction model defined in (|2j. Then 
C(F) = C(F) where C(-) denotes column space. 

Equation |2| shows that the NK algorithm constructs an 
interaction model in an interesting way. Note that the NK 
algorithm dictates that the interaction model contain all sub- 
interactions contained in higher order interactions. For exam- 
ple, if an NK landscape has a fourth order interaction between 
loci {1,2,3,4}, then it also has the four third-order and six 
second-order interactions defined by the subsets of the four 
interacting loci. 

However, knowing which interactions have non-zero coef- 
ficients is not sufficient to fully understand the structure of 
NK landscapes. A complete understanding requires knowing 
the statistical properties of the main effect and interaction 
coefficients. 

A. Induced properties of interaction model coefficients 

While we have established that generalized NK landscapes 
can be represented by interaction models and shown how the 
epistatic interaction sets define the interaction terms included 
in the model, additional insight into landscape properties 
and a complete understanding of the representation requires 
knowledge of the distributional properties of the interaction 
model coefficients. This is addressed in the next proposition, 
with proof given in the appendix. 

Proposition 2. Consider a generalized NK landscape de- 
fined by interaction sets V\ , V2, ■ ■ ■ , V/v and given by p = 
Fw — F/3 with weight vector w where E[wi] = p,/N and 
Var[wi] — a 2 /N. For U G T, let /3jj denote the coefficient of 
the interaction term corresponding to U, and /(•) the indicator 
function, i.e. 

ifU£2 v *; 

1 ifUe2 Vi . 



l(Ue2 v -) = 



Then 



E\ 



E [M =0, for [7^0 
Cov[l3 u ,p u .] = 0, for U^U* 



and 



Var[/3u] 



n N 

-y 

N ^ 

i=l 



2-(^+D /(f/ e 2 ^). 



Some remarks are in order. First, if the weights w are 
normally distributed, then the interaction model coefficients, 
which are linear functions of w, are independent and normally 
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distributed with the indicated means and variances. For other 
distributions on w, the exact distribution of the coefficients 
is an often intractable convolution problem. However, as 
both N and K increase, the central limit theorem can be 
invoked to show that the coefficients will be approximately 
normally distributed (contrast with approximate normality of 
the fitnesses, which only requires N large). Regardless of the 
distribution of w, the coefficients are uncorrelated. 

For w normally distributed, £7[|/3[/|] = yJ\2x[f3ij\2jTT , i.e. 
the expected magnitude of ftu is proportional to it's standard 
deviation. More generally, the variance represents the expected 
magnitude of the square of the coefficient (Var[/3[/] = E[0^j\ 
for U 7^ 0). With this perspective, it is interesting to note 
how the NK algorithm assigns importance (magnitude) to 
the interaction terms. The expected squared magnitude of Pu 
depends on the frequency with which U € 2 Vi . Then for 
the classic NK model where K is constant and K -C N, 
main effects have the largest expected magnitude, second order 
interactions would typically have larger expected magnitude 
than third order interactions and so on. On the other hand, 
when K = N — 1, all coefficients have the same expected 
magnitude because each power set 2 Vi contains main effects 
and interactions of all orders. 

When the coefficients are normally distributed, indepen- 
dence of model coefficients means that, for example, knowing 
the magnitude of a two loci interaction provides no informa- 
tion on the magnitudes of the corresponding main effects of the 
loci-the NK algorithm assigns magnitudes of the coefficients 
completely independently. 

The interaction model representation of a landscape with 
binary loci is equivalent to the representation given by Walsh 
functions [16]. While NK landscapes have been studied from 
the perspective of Walsh functions [6], the above results 
provide a transparent and explicit analysis of the Walsh coeffi- 
cients for generalized NK landscapes, showing exactly which 
Walsh coefficients are nonzero and the statistical properties 
of these coefficients. For example, note that Theorems 4 and 
5 in [6 1 follow immediately from the representation given in 
equation |2]). 

An interesting property of generalized NK landscapes 
gleaned from Propositions [T] and [2] is that two landscapes 
defined by different interaction sets can lead to landscapes 
with the same set of main effects and interactions, i.e. 
identical matrices F, but the variances of the coefficients 
of the landscapes can be different, suggesting that the 
properties of the resultant landscapes may also be differ- 
ent. Consider, for example, the two landscapes defined by 
the following interaction sets: Via — {1, 2, 3, 4}, Vza = 
{2,3},V 3A = {l,3},V iA = {1,3,4},^ = {2,5} and 
V 1B = {IA},V 2B = {1,2,3,4} V 3B = {3},V iB = 
{4} 7 V 5B = {2,5}. Clearly T A = Ujli 2 ^ = U=i 2 VtA = 
T B so that the landscapes have identical column space, but by 
Proposition [2] the variances of the terms are not all identical. 
Consider, for example, the main effect for locus 1: From 
Proposition |2] Var[/3i A ] = 7a 2 /16N > Vai[p 1B ] = 5a 2 /16N. 



B. Rank of NK landscapes and maximal rank models 

Here we show that the rank of the model matrix for an 
NK landscape is determined by the interaction sets {Vi}fL 1 
that define the landscape. We derive the maximum achievable 
rank and describe how interaction sets can be constructed to 
maximize rank. 

Proposition 3. rank(F) = M. (T^ where M(-) denotes 
counting measure and where T = (J i=1 2 *. 

Proof. The result follows immediately from noting that the 
complete set of linear and interaction terms of all orders, 
{l,xi,x 2 , ■ ■ .,x N ,xi*x 2 , ■ • ■ ,x N ^ 1 *x N ,x 1 *x 2 *x 3 , ...,xi* 
x 2 * ■■■ * xn} comprise a Hadamard matrix, a linearly 
independent set of (orthogonal) column vectors. □ 

It follows from the proof of Proposition [3] that rank(F) = 
[number of model interactions — (N + 1)], i.e. model rank is a 
one-to-one function of the total number of interactions in an 
NK landscape. 

The next results give an upper bound on the rank of an NK 
landscape, and provide conditions under which maximal rank 
designs exist for classic NK landscapes. 

Proposition 4. For the classic NK landscape (K constant 
across loci), max{ra«£(F)} < mm{2 N , N2 K+1 + 1 - N(K + 
1)} where the max is over all possible designs for fixed values 
of N and K. 

Proof. The maximum rank occurs when there are no overlaps 
in interactions between the inputs to the loci. In this case 
each locus will contribute £j=V ( K f l ) = 2 K+1 - (K + 2) 
interactions of order two or higher. Summing over the N loci 
and including the intercept and main effects terms gives the 
result. □ 

Proposition [4] begs the question as to when maximal rank 
classic landscapes exist. Equation |2| shows that all classic 
NK landscapes contain a main effect corresponding to each 
loci. Then for both TV and K fixed, maximizing rank is 
equivalent to maximizing the number of interactions in the 
model. From the proof of Proposition [4] it is evident that a 
necessary and sufficient condition is that each set V.- L contain 
unique pairs of loci, as this ensures no redundant two factor 
interactions, and by extension no redundant higher order 
interactions. We now give explicit cases when the maximum 
rank given in Proposition |4] above can be achieved. The proof 
of the following proposition is given in the appendix. 

Proposition 5. There exist classic NK landscapes of maximal 
rank for the following values of N and K. 



1) 


K 


= 2 and all N > 7. 




2) 


K 


= 3 and all N > 13. 




3) 


K 


= 4 and all N > 21 except possibly N — 


22. 


4) 


K 


= 5 and all N > 31 except possibly N = 


32, 33, 34. 


5) 


K 


= 6 and all N > 51. 




6) 


K 


= 7 and N = 57, 64, 67, 69 and all N > 71. 


7) 


K 


= 8 and N = 73, 89 and all N > 91. 




8) 


K 


= 9 and N = 91 and all N > 111. 





It is not difficult to extend Proposition [5] to larger values of 
K, and the approach for doing so is contained in the proof of 
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the result. It is worth noting, and not difficult to prove, that a 
necessary condition for existence of maximal rank landscapes 
is that N > K 2 + K + 1, because otherwise at least one pair 
of loci must occur together in more than one set V*. 

Proposition [5] established the existence of maximum rank 
NK landscapes. We now discuss the construction of these 
landscapes, which is easily achieved through the use of what 
we term NK difference sets. The combinatorial theory under- 
lying these sets is detailed in the proof of Proposition [5] given 
in the appendix. Table [I] provides the needed difference sets 
for up to K — 9. 

To construct a maximal rank classic NK landscape, recall 
that an NK landscape is completely determined by it's in- 
teraction sets {Vi}^ =1 . Interaction sets resulting in maximal 
rank designs are defined by incrementing the elements in the 
difference sets by one (modulo N) until N sets have been 
defined. This process is illustrated in the following example. 

Example: Consider constructing an NK landscape achiev- 
ing maximum rank when N — 7 and K — 2. From Table 
[I] the NK difference set is {0,1,3}. Then define Vi = 
{1,2,4},V 2 = {2,3,5},V 3 = {3,4,6},F 4 = {4,5,7},V 5 = 
{5, 6, 1}, V 6 = {6, 7, 2}, V 7 = {7, 1, 3}. Note that increments 
are modulo N, i.e. increments exceeding N "wrap around". 

TABLE I: NK difference sets used for constructing maximum 
rank designs. The second column gives values of N such that 
maximal NK landscapes exist. 

K N NK Difference Set 

2 N>7 {0,1,3} 

3 iV>13 {0,1,4,6} 

4 iV>23 {0,2,7,8,11} 

5 iV>35 {0,1,4,10,12,17} 

6 iV>51 {0,1,4,10,18,23,25} 

7 iV>71 {0,4,5,17,19,25,28,35} 

8 iV>91 {0,2,10,24,25,29,36,42,45} 

9 TV > 111 {0,1,6,10,23,26,34,41,53,55} 



The design constructed in the example is shown in Fig- 



ure 



la The ith column of the figure gives the elements in V$. 
Notice that no pair of loci appears together more than once 
across columns. The landscape resulting from this design has 
rank N2 K+1 + 1 — N(K+1) = 36, resulting from an intercept, 
7 main effects, 21 two factor interactions and 7 three factor 
interactions. 



Contrast with Figure lb which shows the design resulting 
from choosing adjacent loci for the epistatic interactions. This 
design has rank 29, resulting from an intercept, 7 main effects, 
14 two factor interactions and 7 three factor interactions. 
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(a) Maximal Rank Design 
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(b) Adjacent Loci Design 

Fig. 1: Interaction sets for maximal and adjacent loci designs. 



IV. Number of Local optima 

The number of local optima is perhaps the strongest measure 
of landscape ruggedness and search difficulty [8|, [17|, |19|. 
For fixed N, it has been established empirically that the 
number of local optima increase with K Q. The result 
is not surprising as increasing K increases the landscape 



model rank by increasing both the number and order of 
interactions defining the landscape. An unexplored question 
is the association between the number of local optima and 
landscape rank when N and K are both fixed. 

We begin by providing an analytic expression for the 
expected number of local optima on NK and generalized NK 
landscapes. 
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Proposition 6. Consider a generalized NK landscape de- 
fined by interaction sets V\ , V% , . . . , Vat and weight vector 
w where the elements of w are independent normal random 
variables with E\w^\ = \ijN and Var\u)i\ = a 2 /N. Then the 
expected number of local optima on tlie landscape is given by 
2 JV $(0; E) where 

poo poo poo 

$(0;E)= / / •••/ e-^ T ^ lz d Zl dz 2 ...dz N 
Jo Jo Jo 

and where E is an N x N symmetric matrix with elements 
Oij given by 

{ 2 -wEti^&V k ) ifi=j; 

a ij = l £(Ef=i^eF fe ) + Ef=il(iel4) 

( -TLlHi^Vk or j€V k j) ifi + 3- 

Proof. Let p* denote the fitness of a randomly selected lo- 
cation on the landscape, and let p* denote the fitness one 
hamming distance away obtained by flipping the value of 
locus i. Let z = (p* — p*,p* — pi, . . . ,p* — p* N ) denote the 
N x 1 vector of fitness differences. Then the probability that 
the random location with fitness p* is a local optima is given 
by the orthant probability Pr(z > 0), and the expected number 
of local optima on the landscape is then 2 w Pr(z > 0). When w 
is jointly normally distributed, it follows that z is multivariate 
normal because it is a linear function of w. The mean of z 
is clearly zero, and the variance matrix is straightforward to 
derive and is given by E defined above. The derivation is 
omitted. The orthant probability Pr(z > 0) is therefore given 
by$(0;E). " " □ 

Proposition [6] assumes w is multivariate normal. The utility 
of this assumption is that the expected number of local optima 
then depends on the multivariate normal orthant probability 
$(0;E), allowing us to take advantage of the extensive 
research on the numerical computation of multivariate normal 
probabilities (e.g. and references therein). We are then 
able to estimate the number of local optima for very large 
landscapes without having to generate an actual landscape and 
check individually whether each location is a local peak. Note 
the normality assumption implies that the expected number of 
local optima is a function of only N and E. 

For arbitrary probability distributions for the weights w, 
computation of the orthant probability Pr(z > 0) is typically 
an intractable N dimensional integral. When the weights are 
non-normal and N,K are both large, the distribution of z (see 
proof of Proposition [6]l is approximately normal by the central 
limit theorem, and normal orthant probabilities should then 
provide reasonable approximations to the expected number of 
local optima. 

An empirical study was done to explore the relation between 
landscape rank and the number of local optima for both classic 
and generalized NK landscapes. For each combination of 
N = 25, 50, 100 and K = 1, 2, . . . , 7, we first generated 
interaction sets for 20 classic NK landscapes. For each value 
of N and K, maximal rank landscapes were generated when 
they were known to exist, and an adjacent loci interaction 
design was also generated. Additional classic NK landscapes 



were generated by randomly selecting an N x N latin square 
using the methods in |7| and then randomly selecting K + 1 
rows from the latin square. The resulting N columns of K+ 1 
elements comprise the N interaction sets that define a classic 
NK landscape. For each landscape, the rank and the expected 
number of local optima were computed. Multivariate normal 
orthant probabilities representing the expected number of local 
optima were computed using the "mvtnorm" package in the 
R computing environment, see 0. 

As seen in Figure [2] there is a strong positive correlation 
between the expected number of local optima and landscape 
rank for fixed N and increasing K. Figure [3] provides ad- 
ditional resolution with separate plots for K — 3,4,5 when 
N = 50, showing very strong negative correlation between 
rank and the number of local optima when N and K are both 
fixed. Results for N = 25 and 100 were similar. The maximal 
or near maximal rank designs are given by the points in the 
lower right corners of Figure [3] whereas the adjacent designs 
are given by the point in the upper left corners. Note that 
landscapes computed using adjacent loci to define interaction 
sets resulted in landscapes with significantly lower rank and 
larger expected number of local optima than randomly chosen 
landscapes. 

The strong negative correlation seen in Figure [5] is perhaps 
surprising. Recall that landscape rank is equivalent to the 
number of terms in the interaction model representation, so 
that larger rank corresponds to more interaction terms. It is 
well known that landscapes with main effects but no interac- 
tions have only a single peak. It would seem that additional 
interaction terms would translate on average to more rugged 
landscapes, but this is not the case when N and K are both 
fixed. 

This phenomenon is explained, at least partially, by noting 
that when Proposition [2] is applied to any maximal design for 
an NK landscape, it follows that main effects have variance 
(K + l)2~( K+1 'a 2 /N and interactions of all orders have 
variance 2~( K+1 ^a 2 /N giving a ratio of (K + 1)/1, demon- 
strating that main effects can have considerably more influence 
than interactions in maximal designs. This observation would 
clearly extend to designs that are nearly maximal. Conversely, 
for an adjacent loci design, the variation is spread more 
equitably among main effects and interactions. For example, 
applying Proposition [2] the ratio of variances for a main effect 
and two factor interaction is (K + 1)/K for the adjacent loci 
design. 



Note also from Proposition [2] that for fixed N and K, the 
variance of the main effects and the sum of the variances 
of interaction coefficients is constant for all classic NK 
landscapes, i.e. these quantities do not vary regardless of how 
the interaction sets are defined. It follows that for fixed N 
and K, designs that increase the number of interaction terms 
must have a decreased average variance for the interaction 
terms. In other words, additional model complexity achieved 
through the addition of interaction terms is necessarily offset 
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Fig. 2: Classic 7VK" landscapes: Expected number of local optima versus landscape rank (on log-log scale) for K=\ to 7 



by a reduction in their expected magnitude. 

Additional simulations were done to assess the relation 
between rank and number of local optima for generalized 
NK landscapes. To construct generalized NK landscapes, we 
kept the sizes of the interaction sets Vi constant for a given 
landscape, but did not restrict the number of times that locus 
i could appear in the interaction sets Vj for i ^ j (for classic 
NK landscapes, locus i appears in exactly K interaction sets 
Vi for i + j). 

We randomly generated 20 generalized NK landscapes for 
each combination of TV = 25, 50, 100 and K = 1, 2, . . . , 7. 
A strong positive correlation is again seen between rank and 



expected number of local optima as K varies with N fixed, 
see Figure [4] However, the additional resolution provided by 
Figure [5] shows a positive correlation between landscape rank 
and expected number of local optima. An explanation is that 
for generalized NK landscapes the expected magnitude of 
main and interaction effects can vary, and it is possible to 
have expected magnitudes for a subset of interactions larger 
than those for a subset of main effects. 
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V. Conclusion 

Representation of NK landscapes as linear models in matrix 
form provides a transparent connection to parametric linear 
interaction models, and provides a straightforward means for 
deriving the statistical properties of interaction model coeffi- 
cients induced by the NK landscape algorithm. The interaction 
model representation coupled with distributional properties of 
model coefficients provides new insights into properties of NK 
landscapes. Expressing the expected number of local optima as 
a multivariate normal orthant probability provides additional 
insight into aspects of NK landscapes that affect the number 
of local optima, and also allows for quick computation of the 



expected number of local optima. 

The gaps in the data points seen in the horizontal axes 
of Figure |4] and Figure [5] represent gaps in the ranks of 
NK landscape models and clearly illustrate that classic NK 
landscapes only represent a small subset of possible interaction 
designs, see also for a similar observation. While our 
definition of generalized NK landscapes would fill in some 
of the gaps, the definition still requires that lower order 
interactions that are contained in higher order interactions 
appear in the model. 

The representation of NK landscapes as interaction models 
immediately suggests a more general definition that allows for 
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main effects and interactions of any order to appear together 
without restriction. Such landscapes could be constructed with 
any rank ranging from 1 to 2 , and main effect and interaction 
coefficients could be generated with arbitrary variances and 
correlation structures. In this vein, fl4~l defined a class of 
tunable landscapes allowing general interaction structures and 
studied the performance of GA on these landscapes. 



Using arguments similar to those in Section IV 



it is not 



difficult to show that for a general interaction model with 
normally distributed coefficients, the expected number of local 
optima is represented by a normal orthant probability that 
depends entirely on the variance/covariance matrix of the 



vector of fitness differences. However, it is not obvious how 
to select interactions and assign variances to coefficients that 
will result in landscapes with large numbers of local peaks, or 
more generally that are difficult to search. There are few results 
regarding the size of normal orthant probabilities as a function 
of properties of the variance/covariance matrix, though lfT8ll 
provide potentially useful results in this context. 

Characterizing the search difficulty of landscapes is itself 
difficult |fl~3), and the number of local optima is an imperfect 
determinant of search difficulty [8|, |13|. For example, the 
needle-in-a-haystack landscape, which can be represented as 
a parametric interaction model containing all possible interac- 
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tions where the magnitudes of all main and interactions effects 
are equal, has a single peak which can only be found by 
guessing. While the maximal rank designs for fixed N and 
K had smaller expected number of local optima, they may 
posses other attributes that affect search difficulty. 

An attractive feature of the NK algorithm is that a very 
rich set of landscapes can be generated with only two tuning 
parameters. Whether a procedure with comparable simplicity 
can be developed for constructing a more general class of 
interaction models that are difficult to search is an area for 
future research. 



Appendix 
Proof of Propositions [T]|2]and[5] 

The following is a proof of Proposition [T] 

Proof. We first construct a matrix representation for the inter- 
action model. For the zth locus, let Si t \ (x) denote the vector of 
covariates corresponding to the inputs for locus i. For example, 
if V\ = {1,4,6,9}, then S\,i = (xi, x 4 , x 6 , x 9 ). Note that 
there are Ki + 1 elements in S^i. More generally, let Si )m (x) 
denote the set of m-order interactions for the inputs 

of locus i. For i = 1,...N, define the 1 x 2 Ki+1 vector 
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fi(x) = (l,S it i > S ii 2,--- > S ij ( K+1 )), and note that f ; (-) can 
take 2 Ki+1 different values. 

Consider the 1 x C vector f(x) = {fi(x), f 2 (x), . . .,f]y(x)} 
where C = 2 ^ ,+1 - Define the 2 W x C model matrix 



F* 



( f(xi) \ 
f(x 2 ) 

V f(x 2 w) J 



(3) 



Note that F* =/= F but clearly C(F*) = C(F) as F* 
is comprised of all the columns in F with some columns 
repeated. F* is an overparameterized version of F, and in 
the following we show that C(F*) = C(F). 

The proof uses two basic results. First, suppose two matrices 
have the same column space. Choose I rows from the first 
matrix and form a new matrix by expanding this matrix by row 
concatenating each of the / rows j times. The same operations 
are applied to the second matrix, that is the corresponding 
I rows from the second matrix are repeated j times. The 
expanded matrices are then easily seen to have the same 
column space. The second result is that if C(Ai) = C{B{) for 
I = l,...NthenC(A) = C(B) where A = A\ \A 2 \---\A N 
and B = B\ \ B% \ • • • \ Bn and where | denotes column 
concatenation. 

Let y s = (x^x^,. . -,x iK .) and let yi,i,yi, 2 , • • • ,y i!2 K,+i 
be the binary ordering of yi. Define 



/ 



fi(yi,i) 
fi(yi,a) 



\ 



/ fi(yu) \ 
fi(yi, 2 ) 

\fi(yi,2K+i)/ 



\fi(yi,2K+i)/ 

Then C(a ; ) = C(a £ ) as a t is a 2 K - +1 x 2 K - +1 full rank matrix, 
and a,; is the 2 Ki+1 x 2 Ki+1 identity matrix. Next, let Fj and 

F* denote columns Y?i=\ 2 ^ 1+1 + 1 to Ei=i 2 A ' ;+1 of F and 
F* respectively. Note F t and F* are obtained from a, and a,i 
by repeating rows of these matrices. Then by the first result 
described at the start of the proof, C(F ; ) = C(F?). 

Finally, the result follows from the fact that F* = F* | F£ \ 



Ft and F = F 1 F 2 



F 



N- 



□ 



The following is the proof of Proposition [2] 

Proof. Note that 2 N /? = l T F/3 = l T fw and the expectation 
of the RHS is easily seen to be 2 N fj, from which it follows 
that E[f3q,] = fi. Let F[ k j denote the kth column of F and 
Pu [k] me corresponding interaction model coefficient. Then 
for k > 1, E[f3 U[k] ] = E[Ff k] Ff3] = £[Ff k] Fw] - where 
the last equality follows from l T F[ k ] = whenever k > 1. 

To obtain the variance/covariance results, note that F/3 = 
Fw implies 

P = 2~ n F t Fmv 
where we use the fact that the columns of F are orthogonal, 



each with norm 2 N . Then 

Var[/3] = 2~ 2N F T FF T F x Var[w] 



cr 2 2~ 2N 

2- 2N F T FF T F x Varfwl = F T FF T F. 

N 



To evaluate F T F, let fjj represent the jth column of Fi 
where F{ represents the submatrix of F corresponding to the 
interaction set V, (columns J^lZl 2^ i+1 + 1 to £){ =1 2 K,+1 
of F). Then it is not difficult to show that 



(4) 



where the function h(i,j,k) £ {—1,1} when J7rw E 
2 Vi and zero otherwise (h(i,j,k) is defined explicitly be- 
low), and where Uu.] denotes the interaction term cor- 
responding to F[ k j. Let g k T denote the fcth row of 
F T F. Then it follows from the identity above that 

g k T g k = Ef=i2^ +1 (2 Ar - ( ^ +1) ) 2 ^w e 2^) = 

Eili 2 2JV - (ifi+1) /(J7 [fe ] e 2 Vk ), and then that 

2 N 



1 

N \2~n 



^2 2W -^+ 1 )/(f/ [ „6 2^) 



N 



i=l 

Finally, the covariance result will follow provided the rows 
of F T F are orthogonal. We first define the function h(i, j, k) 
appearing in Q. For j = 1, . . . , 2 Ki+1 , let e^ l {j — 1) denote 
the 1 x 2 Ki+1 vector representation of j — 1 as a binary 
number. The elements of e^ 1 (j — 1), from left to right, 



correspond to values for {xi 1} Xi 2 , . . . , Xi K , 1 }, compare to 
the definition of e^(-) in Section II-B Define h(i,j,k) = 
I(U [k ] G 2^)11^ 2(e- 1 (j - 1) - 1) where the right hand 
side represents the product of the values of the elements of 
e~ x {j — 1) corresponding to the loci in U\\a. For k — 1, we 
have U\k] = Urn = 0, and then define h(i,j,l) — 1 for all 



Example: Consider again the example in Section p-B| where 

N = 3, Ki = 1 for i = 1, 2,3 and V x = {1, 2}, V 2 = {2, 3} 
and V3 = {1, 3}. The interaction model is p — /3@ + fl\X\ + 
P2X2 + P3X3 + /3i2Xix 2 + /3i 3 iiX3 + /3 2 3X 2 x 3 . 

For i = 1, j = 3, e^ x (3 - 1) = (1,0). Then for k = 2, 
f7 [2 ] = {1} and /i(l,3,2) = 1. For fc = 3, C/ [3] = {2} and 
/i(l,3,3) = -1. For fc = 5, C/ [5] = {1,2} and /i(l,3,5) = 
1 x -1 = -1. For fc = 6, E7 [6] = {1,3}, U [6] g 2 Vl and 
h(l,j,6) = for j = 1 to 4. 

The inner product of row fc and ^ of F T F is 

N 2 A 'i + 1 

AT 2 K - + 1 

£ 22J V-2(K 1 + 1) ^ h (iJ lk )h(iJ,l). 
i=l 3 =1 

nKj + l 

Note that Ej=i k)h(i,j, I) = for all i follows from 

orthogonality of Fr k i and Frij. Therefore the rows of F T F 
are orthogonal, and it follows that Cov(/3{/, /9[/» ) = for {/ 7^ 
J7*. □ 
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The following is a proof of Proposition [5] In proving 
Proposition [5] we also define and demonstrate the properties 
of the NK difference sets used to construct maximum rank 
designs. 

We begin with the notion of a packing design or just 
packing for short. These are variants of well-known objects 
from combinatorial design theory, see (T); f° r completeness 
we give the definition here. 

A (N, k)— packing design, consists of a set S on N elements 
(called "points") and a collection of subsets of S (called the 
"blocks") all of size n = K + 1 with the property that 
each pair of points in S is contained in at most one of the 
blocks. To obtain an NK landscape of maximum rank, one 
can use the points and the blocks of a (N, k) packing design 
to construct the N interaction sets defining the landscape. The 
next theorem gives this connection. 

Proposition 7. If there exists an [N, k) packing design, then 
there exists a classic NK landscape with N loci and K = 
K — 1 which achieves maximum rank. 

Proof. As noted in the proof of Proposition |4] maximal rank 
designs will occur when each pair of loci occurs at most once. 
The existence of a (N, n — K + 1) packing design ensures 
that there are no redundant two factor interactions in the NK 
interaction set specification, and by extension no redundant 
interactions of order greater than two. Hence a NK landscape 
of maximal possible rank will result. □ 

In order to construct NK landscapes via Proposition [7] we 
must construct (N, k) packing designs. Constructing (TV, k) 
packing designs is made easier by employing so called differ- 
ence methods from combinatorial design theory. Let Z n be the 
integers modulo n, so Z„ = {0, 1, 2, ... n — 1} with addition 
modulo n. Define an (N, k) difference set in Zn to be a subset 
D = {xx, X2, ■ ■ ■ , Xk,} C Zjv of size k with the property that 
the list of all differences Xi — Xj (where x.- L , Xj £ D and i ^ j) 
contains each nonzero element at most one time. 

An example of a (8, 3) difference set is D — {0, 1, 3}. Note 
that in Z s , 1 - = 1, 0-1 = 7, 3-0 = 3, - 3 = 
5, 3—1 = 2, and 1 — 3 = 6 and hence all of the differences 
are different as required. 

We can use a (TV, k) difference set in Zn to construct an 
(N, k) packing design and thus via Proposition [7] we will have 
constructed a classic NK landscape with N loci and K = k—1 
which achieves maximum rank. The construction is given in 
the next proposition. 

Proposition 8. If there exists an (N, k)— difference set in Zn, 

then there exists an (N, k) packing design. 

Proof. Let D = {x\, X2, ■ ■ ■ , x K } be an (iV, k) difference set 
in Z N - For g e Z N define D+g — {xi+g, x 2 +g, ■ ■ ■ , x K +g}. 
D + g is called a translate of D. We claim the set of all 
translates of D (namely {D + g : g G Zn}) form the blocks 
of an (N, k) packing design. First note that there are exactly 
N of these translates (one for each element in Zn) and also 
note that there are N points (the elements of Zn). We must 
show that any pair of elements a,b e Zn that a and b are in at 
most one of the translates of D. Assume to the contrary that the 



pair {a, b} occurs in two of the translates, say {a, 6} C D +p 
and {a, b} C D + q where p ^ q are both elements of Zm- 
Then a = Xj + p and b = Xj + p for some i and j and 
also a — Xk + q and b = x m + q for some fc and m. Hence 
a — b = Xi — Xj — Xk — x m . From the difference property 
of D this implies that i = k and j — m which in turn says 
that p = q, a contradiction. Hence we have shown that the n 
translates of an (N, n) difference set in Z^ give the blocks of 
an (TV, k) packing design. □ 
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